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Abstract 



In this paper we analyze the convergence of the splitting method for 
shallow water equations. In particular we give an analytical estimation of 
the time step which is necessary for the convergence and then we study the 
behaviour of the motion of the shallow water in the Venice lagoon by using the 
splitting method with a finite element space discretization. The numerical 
calculations show that the splitting method is convergent if the time step of 
the first part is sufficiently small and that it gives a good agreement with the 
experimental data. 

1 Introduction 

The equations of a Newtonian (viscous) fluid are called Navier-Stokes equa- 
tions and in the case of a inviscid fluid they give the Euler equations. The 
Navier-Stokes equations are given by 

+ (u ■ V)u + Vp - z/V 2 u = f , 

at 

V-u = 0, (1) 
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where p = P/p is the reduced pressure and v = p/p is the reduced or 
kinematic viscosity 1 . 

The Navier-Stokes equations are to be integrated in the space-time do- 
main f2x]0,T[c R 3 x R + once an appropriate set of initial and boundary 
conditions has been defined. While the former conditions, in general, only 
need to reproduce a feasible shape of the current solution, the latter are to 
be set and treated with much care due to the effect they can have on the 
evolution process 2,3 . 

In the shallow water hypothesis one assumes that the characteristic hor- 
izontal scale L for the motion (the wave length) is longer than the average 
height h of the fluid, i.e. L » h. In such hypothesis vertical acceleration 
and velocity are negligible and the flux becomes almost horizontal 4 . 

Definition 1.1 The total height of the water is given by 



where H is the depth of the water in the stationary condition above the ref- 
erence level x 3 = 0, and i] is the time-dependent difference, i.e. the height of 
the free surface. 

The external forces acting on the water are of extreme importance to 
determine the equations of motion of the system. 

Definition 1.2 The force of gravity is defined as 



where g = [0, 0, — g] T is the vector acceleration of gravity. 

If we consider large regions of water it is necessary to include other forces, 
like the Coriolis force and the Chezy force, which models the friction of the 
water at the bottom. 

Definition 1.3 The Coriolis force is defined as 



where uj = [0, 0, W3] T is the rotation vector of the earth and = k is called 
Coriolis coefficient. 



h(xi, x 2 , t) = H(xi, x 2 ) + r)(x 1 ,x 2 , t) , 




■cor 
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Definition 1.4 The Chezy force is defined as 

p f 9u\u\ 



where \u\ = \Ju\ + u\, g is the scalar acceleration of gravity, h is the total 
height of the water and k\ is the Chezy coefficient. 

The first step to obtain the shallow water equations is to put u 3 and du 3 /dt 
equal to zero in the Navier Stokes equations. Then the third equation can 
be integrated between —H and rj. The system is yet 3-dimensional because 
the velocities U\ and u 2 are functions also of the x 3 variable. It is possible to 
obtain a 2-dimensional system by performing the substitution 

u 1 (x 1 ,x 2 ,x 3 ,t) -> a(tp)u 1 (x 1 ,x 2 ,t) , 

u 2 (x 1 ,x 2 ,x 3 ,t) -> a(ip)u 2 (x 1 ,x 2 ,t) , 



where 



ip = — with / a(ip)dip 

" Jo 



r] + H 

and then by integrating the new equations over the x 3 variable. In this way 
one can prove the following theorem (see Ref . 5 and 6 for details) . 

Theorem 1.5 Let us consider the inviscid shallow water with external forces 
given by the gravity force and the Coriolis force. The 2-dimensional viscid 
shallow water equations are 

dt dxi dx 2 

where U = [rj, ui, u 2 ] T is the vector of the conservative variables, the flow 
vectors F 1 and F 2 are 

F 1 = [Hu u g-n, 0] T , F 2 = [Hu 2 , 0, grf , (3) 

and the source vector R s is given by 

T5 rn ; 9Ui\u\ gu 2 \u\ T 

R s = [0, k u 2 - -grff, -koui - -j^jf] , ( 4 ) 

where g is the acceleration of gravity, ko and k\ are the coefficients of Coriolis 
and Chezy and H(xi, x 2 ) is the depth of the water in the stationary condition. 
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The shallow water equations are to be integrated in the space-time do- 
main f2x]0,T[c R 2 x R + once an appropriate set of initial and boundary 
conditions has been defined 

U(xi,x 2 ,0) = U (xi,x 2 ) V(x 1 ,x 2 )gO, (5) 

U(x 1 ,x 2 ,t) = V ga (x 1 ,x 2 ,t) V(x 1 ,x 2 )EdSl,Vte]0,T[. (6) 

While the initial conditions, in general, only need to reproduce a feasible 
shape of the current solution, the boundary conditions are to be set and 
treated with much care due to the effect they can have on the evolution 
process. 



2 The splitting method 

A semi-implicit method is adopted for the solution of the system (Q), based 
on the following splitting (see Ref. 7 and 8) 

Fi = F* + F** ,2 = 1,2 (7) 

and 

U(n+1) = u(n) + AU * + AU ** ; (g) 

where AU* and AU** are the increments of the solution vector. 

In the iterative scheme we put F* = so that Fi = F**, i = 1, 2, and the 
system (|j) can be divided in: 

<9AU* 



dt 

and 



R s ■ (9) 



<9AU** dF** dF** , . 

— + a^ + l^ = - (10) 

These equations are integrated in time, in turn, by using an explicit Taylor- 
Galerkin method 8 for (Q) and an implicit ^-method 9 for (|I0|). 

The equation (^|) is discretized in time by using a Taylor expansion to the 
second order 
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where r is the time step. From 

dAU* _ d 2 AU* _ dR s dR s dAV* 

~^r~ Ks and ~d^~~df~dAu* dt ' (12) 

the equation fll!]) can be written 

(AU*) (n+1) = r(R s )^ + y (GR s ) (n) , (13) 

where 

e = <^- <"> 

Because of the computational complexity in the evaluation of the right term 
of equation ( JEf ) , we use a two-step version of the Taylor-Galerkin algorithm. 
This is given by an approximation of (\jy n+1 / 2 ) and R s *- ra+1//2 - ) the Taylor 
expansion at step {n + 1/2) 

U(n+l/2) = u(n) + l( Rs )(«) f (15) 

(R s )(" +1 /2) = (Rs) W + r fdB^ (n) = (Rs)(n) + ^ (GRs) (n) ; (16) 

from which we obtain (GR S )^ as 

(GR s ) (n) = - [(R s ) (n+1/2) - (R s ) (n) l . (17) 

It follows that (AU*) (rt+1) can be written as 

( AU *)(n+l) = r ( Rs )(«+i/2) ; ( lg ) 

and the explicit scheme of the splitting method results 

U(n+1) = jj(n) + R (n+l/2) _ (lg) 

The equation (|T0| ) is given, in explicit form, by 

I (At?") + afr i H *) + A = 

I ( Am D + at W = (2°) 

I (AO + & W = • 
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The time discretization is obtained by using the ^-method 



(Ar/**) (n+1) 
(A<) (n+1) 



+ r 








(21) 



where r is the time step, #i and #2 are real parameters between and 1 and 



u 



(n+fli) 



*\(n+l) 



+ (Am 



**\ (n+l) 



^(n+e 2 ) = ^(n) _|_ 5) 2 (A^**)(" +1 ) 



1,2 



(22) 



Note that the term (Ar;*)( n+1 ) does not appear because the splitting method 
does not include variations in the water elevation due to fll3|) . Substituting 
(52) into pT|), we can write 



(A?T) (n+1) +^iELi^ # (A«; 



,**\(n+l) 



(A«; 



**\(n+l) 



(A M **) (n+1) 



(Ar/**) (n+1) 



(At; 



**\(n+l) 



-T9£- 2 v {n) 



(23) 



Here (Aw**) (n+1) and (Au**) {n+1) can be obtained from the second and third 
equation as a function of (Ar]**y n+1 ' and then substituted in the first equa- 
tion which becomes 



(At; 



~2 



(A77* 



— r 



&L 1 £:[H(uM + e 1 (Au*t 



+1) 



r^ELi£"(^ (n) )} • 
(24) 



3 Time Convergence of the splitting method 

In this section we study analytically the time convergence of the the split- 
ting method. Our study is motivated by the stability analysis of multilevel 
methods for the numerical simulation of turbulent flows performed by Roger 
Temam. He analyzed a simple model, namely a pair of ordinary coupled dif- 
ferential equations, by using numerical schemes based on different treatments 
and different time steps for the two variables of the system 10 . 
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In our splitting method, starting from the time discretization of the pre- 
vious section, we have 



v (n+l) = ^(n) + (ATf + (Ar]**) (n+1) 

v!( l) + (AulY n+1) + (A<) (n+1) 



(25) 



u ( 2 n) + (Au* 2 Y n+1) + (Au* 2 *Y n+1) 



We observe that the first component of R^ n+1 / 2 - ) is zero so that (Ar]*)^ n+1 ^ 
It follows that the first part of the splitting scheme reads 



^(n+l) = ^(n) 



u (n+l) = u (n) + (Am * 



*\(n+l) 



where 



(Aul)^ = au^ + (3u { 2 n) 



(26) 

(27) 
(28) 



with a = (-^tt - r D + ^), /3 = (rA; - r 2 /J) and /J = #|u|/(£; 2 # ). Here 
we suppose that D is constant during the iterative process. The iteration 
scheme can be written as 



ut +1) 



( \ 



j* 



(n) 

(n) 
«2 



(29) 



where 



J* 



/ 1 

1 + a (3 
\ -/3 1 + a / 



(30) 



Remark 3.1 We observe that the first part of the splitting scheme and its 
matrix J* do not depend on Q\ and 9 2 . 

Let us study the properties of the first part of the splitting scheme gen- 
erated by the matrix J*. 



7 



Definition 3.2 The spectral radius of the matrix J* is given by 



P(J*) 



max{\\i\, i = 1,2, 3} 



where A« 7 i — 1,2, 3, are the eigenvalues of the matrix J* . 

By solving the eigenvalue equation we find that the eigenvalues of J* are 
Ai = and A 2 ,3 = (1 + a) ± i(3. 

Definition 3.3 The iteration scheme induced by the matrix J* is called con- 
vergent if the spectral radius is such that p(J*) < 1. 

Now we can prove the following theorem. 

Theorem 3.4 The first part of the splitting method, induced by the matrix 
J* , is not convergent for D = 0. For D ^ it is convergent if and only if 
the following inequality is satisfied 



{k* + D 2 (A - kl) + 4D 2 )r 3 - 4-D(-D 2 + 2k - k 2 )r 2 + 8D 2 t - 8D < , (31) 



where D = g\u\/ (k\H). 

Proof. The convergence condition for the iterative process is such that the 
spectral radius of J* is less than one. We have seen that the eigenvalues of J* 
are Ai = and A 2j 3 = (1 + a) ± if3. It follows that the convergence condition 

is (l + a) 2 + (3 2 < 1, where a = (r 2 f -tD + t 2 ^) and f3 = (Tk -r 2 D), and 
this condition gives the inequality of the theorem. In particular, if D = we 
have the inequality /cq < which is never satisfied. □ 

Remark 3.5 The convergence of the iteration scheme generated by the ma- 
trix J* does not depend on f, 9\ and 62 because these parameters do not 
appear in the eigenvalues of J* and in the inequality of the theorem 3. 3. 

The equation associated to the inequality of the theorem 3.3 is given by 



ar 3 - br 2 + cr - d = , 



(32) 
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where a = (Jfeg + D 2 {A - fcg) + 4L> 2 ), 6 = 4D(£> 2 + 2A; - fcjj), c = 8L> 2 
and d = 8D. This equation has two complex conjugate solutions and a real 
solution given by 

r - A _ 2l/3 (~ fe2 + 3ac ) + * V3 (33) 
Tc "3a 3agV3 + 31/33^ > ^ 



where q = 2b 3 - 9abc + 27a 2 d + ^4(-6 2 + 3ac) 3 + (26 3 - 9abc + 27a 2 d) 2 . It 
follows that the iteration scheme induced by the matrix J* is convergent if 
and only if the following inequality is satisfied 

r<r c . (34) 

In particular, with the positions g = 9.81 m/sec 2 , k = 1CT 4 sec^ 1 , k\ = 40 
m 1 / 2 sec~ 1 , \u\ = 0.1 m/sec and H = 0.1 m, we find r < t c = 5.41 sec. 
We shall use this estimation for the numerical implementation of the splitting 
algorithm. 

In conclusion, we see that the convergence of the first part of the split- 
ting method depends strongly on the time step r. This time step must be 
very small to ensure the convergence of the iteration scheme. Moreover the 
inclusion of the Chezy force is essential for the convergence. In fact, if the 
Chezy force is not included (D = 0) then the iterative process diverges. 

We observe that the problem of spatial discretization can be studied 
with the Finite Element Method (see next section). As shown by Ciarlet 
and Raviart 11 , the interpolating function of Finite Element Method controls 
completely the spatial convergence. It follows that, to analyze the time con- 
vergence of the second part of the splitting method, the spatial dependence 
can be neglected. In this way, we get (Am**)^ +1 ^ = 0, % — 1,2, and the 
splitting scheme reads 

jjin+l) _ ^(n) 

u (n+l) = u {n) + ( Au *)(n+1) (35) 

u< n+1 > = U W + (A«;)< n+1 > 

where 

oxi oxi ox 2 oxi ox 2 0x2 
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*\(n+l) 



(n) a (n) 



(A M ;)(« +1 ) = -/?4 n) + «4 n) 



(36) 



with a 



-Zp. -rD + ^f-), {3 = (rk - r 2 D) and D = g\u\/(k\H). It 
follows that the iteration scheme can be written as 



I v (n+l) \ 



ut +1) 



J** 



( \ 
) 



(n) 
in) 



(37) 



where 

J** : 



1 + a 



/3 
1 + a 



(38) 



Remark 3.6 TTie iteration scheme induced by the matrix J** does not de- 
pend on 02- 

It is easy to see that the eigenvalues of J** are the same of J*. They 
are Ai = and A2,3 = (1 + at) ± i/3. It means that also the convergence of 
the iteration scheme generated by the matrix J** does not depend on f, 9\ 
and 82 because these parameters do not not appear in the eigenvalues of J**. 
In conclusion, we do not have limitations for f. In practice, because of the 
effect of the error propagation due to the numerical approximations, we use 
f ~ 5 minutes. 



4 Finite element space discretization 



The spatial approximation with the finite element method is obtained by 
using linear form functions for the integer step (n, n + 1, . . .) and constant 
functions for the half step (n — 1/2, n + 1/2, . . .). In this way the equation 
reads (see also Ref. 12 and 13) 



(MAU*)} 



(n+l) 



T 



(R s ) (n VA4l - (R s ) (n+1/2) - (RJ (n) 
z Jn r 1 



{(R-) 



(n+l/2) 



+ 



;R s ) (n) -(R s ) (n) ]} <Mx 



j 



1,2, 



bjdx. (39) 



10 



where N is the total number of nodes, <j>j is the weight linear function of the 
node j, the bar denotes mean values calculated on the element and 



M 



[M]ij = J Mjdx i, j = 1, 2, ■ ■ ■ , N (40) 



is the mass matrix. 

By using again linear triangular elements for the spatial discretization 
with the finite elements method, and by using the Green formulas to the 
terms which include the second derivatives, the variational formulation (f24|) 
gives the following system of discrete equations 

(M + fVi^S)(Ap**)(™ +1 ) = 

- f {eLi Q* [h (*4 n) + 0! (A M *) (n+1) )] + fOiSp^} 



where 



i=l 



When (Ar/**)^ 1 ^ is evaluated by fl4l|), the second and third equations of 
(gH) can be solved in (Awf) (n+1) and (Au* 2 *) {n+1) '; their discretization to 
finite elements is given by 

M(A<*) (n+1) = -fgQi U n) + 9 2 (Ari**) {n+1) ] % = 1, 2 . (43) 



It is easy to see that this set of equations is similar to Q3"9p; as a consequence 
it can be solved in the same way. 

In summary, the solution at each time step implies the following opera- 
tions: 

a) solve (H) for (Au*) (n+1) and (Au*) (n+1) (here (Arf )( n+1 ) = because 
the first component of R s is zero); 

b) calculate (Arf*)(" +1 ) from (0); 

c) solve Qg) for (Aw**) (n+1) and (Aw**)( n+1 ); 

d) calculate Tj( n+1 ) = U'"' + (AU*)( n+1 ) + (AU**)( n+1 ). 
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5 Numerical calculations 



In this section we describe our calculations of the motion of the water into 
the Venice lagoon 12 ' 13 . The equation considered are the shallow water equa- 
tions. All simplification are obtained considering different depth of the lagoon 
and in particular a finer discretization is considered where the water runs 
with greater speed namely along the more important canals. The numerical 
method is based on the splitting scheme and the finite element discretization 
presented in the previous sections. 

To analyze the Venice Lagoon we must add the wind effect. Let v the 
wind velocity respect the water and £ an adimensional constant, the wind 
term is given by £\v\v/H. In the approximation that the water velocity is 
negligible respect the wind velocity, the wind term is an external parameter to 
the system, i.e. the wind velocity respect the soil. We take a time-dependent 
wind velocity but space-independent 14 . 

Following our analytical estimation of the time step which is necessary 
for the convergence of the splitting scheme, we choose a time step r = 3 
seconds for the explicit part and f = 5 minutes for the implicit part. In this 
way the model is very stable and capable to preserve the water quantity of 
the system. We decompose the Venice lagoon into 10 subdomains and the 
nodes of the interfaces collected into a unique set. The number of nodes 
of the entire discretization is 1967 and the elements are 3423. The results 
of the calculation obtained using the domain decomposition method without 
parallelization have approximately a speed-up of 25 per cent over a Conjugate 
Gradient solver used on a not-decomposed domain 12,13,17 . 

The model has been verified by using the experimental data of the Uflicio 
Idrografico and Mareografico di Venezia (month of Fabruary 1994) 15 . The 
data have been introduced at the mouth of the lagoon as boundary conditions 
on the open boundary. Therefore a comparison has been made between the 
simulation and the real data on the internal measurement stations of the 
lagoon 16 . The experimental and simulated data present a good agreement in 
all the internal measurement stations 12,13,17 . 
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